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Thanks to their high quality factor, combined to the smallest modal volume, defect-cavities in photonic 
crystal slabs represent a promising, versatile tool for fundamental studies and applications in photonics. In 
paricular, the L3, HO, and HI defects are the most popular and widespread cavity designs, due to their 
compactness, simplicity, and small mode volume. For these cavities, the current best optimal designs still 
result in Q- values of a few times 10 5 only, namely one order of magnitude below the bound set by fabrication 
imperfections and material absorption in silicon. Here, we use a genetic algorithm to find a global maximum 
of the quality factor of these designs, by varying the positions of few neighbouring holes. We consistently 
find Q-values above one million - one order of magnitude higher than previous designs. Furthermore, we 
study the effect of disorder on the optimal designs and conclude that a similar improvement is also expected 
experimentally in state-of-the-art systems. 

Photonic crystal (PhC) cavities are a quintessential element of integrated photonic devices. Thanks to ultra- 
high quality factors and mode volumes close to the diffraction limit 1 " 13 , these devices hold promise in a wide 
range of applications including non-classical light generation, all-optical computational paradigms, solid- 
state cavity quantum electrodynamics, and sensing 14 24 . Starting from a two-dimensional PhC consisting of a 
lattice of air holes etched in a dielectric slab, an optical cavity can be created by introducing a point-like defect - 
e.g. one or few missing or shifted holes - in the otherwise periodic structure. 

A major effort during the last decade has been devoted to the optimization of these structures, in particular 
through the maximization of the quality factor Q and the minimization of the volume V of the cavity mode, as 
optical nonlinearities, Purcell effect, and radiation-matter coupling all depend directly on Q and inversely on 
V 25 28 . To this purpose, three different approaches can be broadly defined. The first is the inverse problem 
approach 3,4 , where an effective equation for the dielectric profile is defined starting from the desired shape of 
the cavity mode, through a semi-analitical formalism. The second is the topology optimization method 29 , where 
variations of the entire topology of the PhC are allowed, and the objective function (either QorQ/V) is maximized 
numerically 30 32 . Although some of these works 3,32 have achieved remarkably high values of Q and Q/V, the 
resulting cavity designs often pose serious technological challenge in terms of manufacturability, as they present 
excessively small holes or holes with irregular pattern and sharp features. A third and completely different strategy 
consists in optimizing simple PhC cavity designs by tweaking only a few geometrical parameters (e.g. by shifting 
the positions and varying the radii of nearby holes), in order to preserve the small spatial footprint and ease of 
fabrication of the design. This approach has produced encouraging results, in particular for the three most 
widespread cavity designs, namely the L3 cavity 1,2,9,33 , the HO cavity (also known as "zero-cell" or "point-shift" 
cavity) 6,9 , and the HI cavity 8,9,11,13,34 . It brought in some cases an increase of the quality factor by more than one 
order of magnitude - reaching values of a few hundred- thousands, or even above one million for the H 1 hexapole 
mode 11,13 - while the mode volume was only slightly increased or sometimes even reduced. A common feature of 
all these optimization works however is the lack of an exhaustive exploration of the parameter space in search of a 
global maximum of the objective function. 

Experimentally, the quality factor of PhC cavities is limited by extrinsic losses, due to absorption and fabrica- 
tion imperfections > 2 - 35 - 37 . More precisely, the measured quality factor Q e can be expressed as 

1 _ 1 

Q e ~ Qt 



l l 



(i) 



where Q t is the theoretical quality factor expected from the ideal structure, while 1/Q„ and lIQd are measures of 
the additional loss rates due to material absorption and to disorder-induced extrinsic losses respectively. For 
silicon PhCs and wavelengths in the 1.5 fim range, record values of Q e ranging between one and five million were 
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measured on cavity designs based on a modulation of a one-dimen- 
sional PhC waveguide 5-7,10-12 . For these designs, Q t ranges between 2 
X 10 7 and 10 s , suggesting a value of several millions for both Q d and 
Q a . While displaying a mode volume up to three times larger than 
that of defect cavities, these waveguide-based designs have a consid- 
erably larger footprint. There are no fundamental reasons that 
should prevent the quality factor of the most popular defect cavities 
from reaching values of several million, close to the current bound set 
by disorder and absorption. 

Here, we adopt a simple optimization strategy. Similarly to several 
existing workg 1 - 2 *** 11 - 13,3 *^ we choose a small set of variational para- 
meters - typically the spatial shifts of a few holes next to the defect - 
thus producing designs that can be easily realized with current nano- 
fabrication processes. Differently from these works however, we 
carry out a global exploration of the parameter space by means of 
an evolutionary algorithm. Here, by "global" we mean the exhaustive 
search for the global maximum of the quality factor in the parameter 
space of choice. In this way we demonstrate that it is possible to 
systematically optimize L3, HO, and HI cavities to values of Q t well 
above 10 6 - typically more than one order of magnitude above pre- 
vious optimal values - without a large increase of the mode volumes. 
The key to this drastic improvement is the exhaustive search, that 
finds configurations overlooked by previous approaches, as exem- 
plified by the simplest of the two L3 designs considered here 2 . This 
optimization procedure is made computationally feasible thanks to 
the use of the guided-mode expansion (GME) method 38 , that allows 
calculation of the modes and quality factors of each variation within 
minutes of computational time. Thanks to this computational 
advantage, we also statistically analyze the influence of fabrication 
imperfections on the optimal designs and conclude that a consid- 
erable improvement in the experimental quality factor can be 
expected, as recently demonstrated for the HO design 39 . 

Results 

The cavities studied here are formed in a triangular lattice with pitch 
a of air-holes of radius R etched in a silicon (n = 3.46) slab of 
thickness d. In what follows, all lengths will be expressed in units 
of a, as the quality factor is invariant upon a spatial rescaling of the 
structure. We however set parameters such that, for the typically used 
thickness d = 220 nm, the resonant modes lie in the telecommuni- 
cation band around A = 1.55 ^m. For the simulation of a single 
structure, we use the GME method 38 (see Methods). This method 
has already proven reliable for modelling high-Q cavities 35,37 . As a 
further check of its validity, the results for all final (optimized) struc- 
tures are verified using the 3-D finite-difference time-domain 
(FDTD) method 40 (see Methods). For each of the optimized designs 
presented here, we also analyse the probability density of Q e in the 
presence of fabrication imperfections (and neglecting the absorption 
loss contribution 1/Q„). The disorder model is that of Gaussian dis- 
tributed random fluctuations in the position and radius of each hole, 
with zero mean and standard deviation a (which is a measure of the 
disorder magnitude). Disorder in all holes in the computational cell 
was included, and no hole-hole correlations were taken into account. 
Irregular hole shapes can in general easily be included, but the result- 
ing effect is well described by an effective radius fluctuation 41 . The 
disorder model thus captures well what is accepted as the main 
source of losses in silicon photonic crystal cavities 12,35,36 . 

The first cavity design we investigate is the widely- employed L3 
cavity formed by three missing holes in a row (Fig. 1(a)), with d = 
0.55a and _R = 0.25a. The quality factor of this cavity has already been 
optimized 2 with respect to shifts of the positions of three neighbour- 
ing holes (marked S ix , S 2x , and S 3x in the Figure), by using a simpli- 
fied approach in which, starting from the unshifted design, each of 
the three shifts has been varied once while keeping the two others 
constant. To explore the extent to which this approach is suitable, we 
compute a full map of the quality factor on a relevant region of the 



(Six, S 2x , SsJ-space. The map is displayed in Fig. l(b)-(m). There, in 
each panel, Q t is plotted as a function of S 2x and S 3x , while the value of 
S ix increases from 0.15a to 0.37a in steps of 0.02a when going from 
panel (b) to panel (m). Technically, these plots already provide a 
global optimization of the cavity (a clear maximum of Q t can be 
identified), although performed in the least practical, brute-force 
way. If applied to the panels of Fig. 1 (though approximately, given 
the coarse Si x step used in this figure), the simplified optimization 
procedure in Ref. 2 leads to the point marked by a white cross in panel 
(e) (more precisely, S lx = 0.21a, S 2x = 0.01a, S 3x = 0.23a), i.e. far off 
the maximum that can be seen in panel (k) at S ix = 0.33a, S 2x = 
0.26a, S 3x = 0.10a. It is also interesting to note that within the range 
of the plots in Fig. l(b-m), two maxima are visible - a local one 
around S lx = 0.25a, S 2x = 0.09a, S 3x = 0.21a in panel (g) and another 
one which first appears at S lx = 0.29a, S 2x = 0.18a, S 3x = —0.10a in 
panel (i) and then shifts to become the global maximum at S ix = 
0.327a,S 2x = 0.257a, S 3x = 0.116a. In general, an even larger number 
of local extrema might in principle be present, especially for larger 
number of parameters, thus making the search for the global 
extremum more difficult. This highlights the need for a global optim- 
ization procedure instead of a more conventional algorithm (e.g. the 
conjugate gradient) that would almost inevitably find a local rather 
than a global maximum. 

Ideally, we would like to apply a global, stochastic (since there is no 
general way to come up with a "good" guess for a starting point) 
procedure to the problem of optimizing the cavity parameters. Thus, 
we choose to employ a genetic algorithm which typically relies on a 
range of evolution-inspired techniques to create consecutive "gen- 
erations" of cavities each containing better and better "individuals", 
until convergence is reached. This type of algorithm proves very 
efficient especially with increasing number of free parameters. We 
use the genetic algorithm provided in the MATLAB® Global 
Optimization Toolbox (see Methods), where we choose the objective 
function to be the GME-computed quality factor Q t . When applied to 
the L3 with freedom in S lx , S 2x and S 3x , the optimal design is found 
for Si* = 0.327a, S 2x = 0.257a and S 3x = 0.116a (Fig. 2(a)). This 
yields Q t = 2.1 X 10 6 (FDTD: 1.6 X 10 s ), which is an increase by a 
factor of ~ 6 as compared to the previously highest value 2 of 3.3 X 
10 5 (FDTD: 2.6 X 10 5 ), while the mode volume (see Methods) 
increases from 0.77(A/«) 3 to 0.94(A/«) 3 . One obvious advantage of 
using evolutionary optimization rather than the brute-force para- 
meter scan of Fig. 1 is the precision with which the maximum can 
be pinpointed; another one is that a few tens of generations with a 
population of 80 individuals are sufficient to reach this converged 
configuration (see Methods). Moreover, the design can be further 
improved if two more shifts (S 4x and S 5x in Fig. 1(a)) are allowed in 
the optimization, which is still easily handled by the genetic algo- 
rithm, although —100 generations are needed for convergence. In 
this case, the optimized design is found for Si x = 0.337a, S 2x = 
0.270a, S 3x = 0.088a, S 4x = 0.323a, and S Sx = 0.173a (Fig. 2(b)), 
and has Q t = 5.1 X 10 6 (FDTD: 4.2 X 10 6 ) with mode volume V = 
0.95(/L/n) 3 , i.e. an increase in Q t by one order of magnitude compared 
to the previous optimal values 2,33 , with an increase in the mode 
volume comparable to the three-shift case. The resonant frequency 

of the modes is at =0.259 for both designs, which is slightly 

2ttc 

lower than the frequency- — = 0.263 of the unmodified L3 (i.e. with 

the same dla and Rla but with no hole shifts). 

The choice of just these few free parameters is not unique, but 
appears optimal. Attempts were made at including other position 
shifts or radii variations of the holes surrounding the cavity defect, 
which however did not bring to a significant further increase in Q t . 
Having only a few hole shifts as free parameters results in a tech- 
nologically friendly structure, and in a more compact cavity defect, 
characterized by a much smaller footprint on the PhC, than wave- 
guide-based ultrahigh- Q designs 4,5,7 . In addition, the present designs 
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Figure 1 | (a): The design of the L3 cavity. For quality factor optimization, shifts of the positions of the five neighbouring holes in the x-direction 
were introduced, labeled as Si„ S 2x , S 3x) S 4x , and S 5x in the figure, (b)-(m): A parameter scan of the GME-computed quality factor values for different 
Si» S 2x and S 3M where S ix starts from 0.1 5a in panel (b) and increases in multiples of 0.02a in every consecutive panel, up to 0.37a in (m), and S 4x = S 5l = 0 
in all panels. 



are as robust to fabrication imperfections as any other high-Q PhC 
cavity 37 , as can be inferred from Fig. 2(c)-(f). In panel (c) we plot, for 
the three-shift L3, the dependence of Q t on S 2x and S 3x as in Fig. 1, but 
for the value S lx = 0.327a corresponding to our optimal design (the 
white cross indicates where the design lies with respect to S 2x and 
S 3x ). In this plot, we observe that the width of the maximum is larger 
than the typical uncertainty in the hole positions (smaller than 
0.005a for Si 42 ). Furthermore, in Fig. 2(d)-(e) we show, respectively 
for the three- and five-shift design, the computed probability of 
occurrence of Q e - values in presence of disorder. Each of these histo- 
grams was obtained by simulating 1000 disordered realisations of the 
corresponding cavity design. The blue plot in panel (d) in particular 
shows that for a state-of-the-art disorder magnitude a = 0.0015a (i.e. 
about 0.6 nm 12,35 , when assuming a = 400 nm in a silicon slab), the 
average value lies at about Q e = 2.5 X 10 6 , i.e. quality factors one 
order of magnitude larger than the previous theoretical maximum 
can be expected in practice, highlighting the significance of the design 



optimization. Finally we note that, for a given set of optimal values of 
the S„ x parameters, the designs are also robust to small changes in the 
overall hole radius _R and slab thickness d, that can originate from an 
offset in the fabrication process and/or be introduced on purpose in 
order to e.g. tune the resonant frequency to a desired value. To show 
this, in panel (f) we plot the value of Q t obtained by varying R and d 
while keeping the shifts S lx — S 5x constant and set to the values 
obtained for the optimal design computed at d = 0.55a and R = 
0.25a. We observe that Q t > 4 X 10 6 for a range of _R and d values 
which is much larger than the fabrication uncertainty and which 
allows fine-tuning of the frequency. For certain values of d and R 
(to the right of the dashed line in the Figure), higher-order guided 
modes of the slab become non-negligible. More precisely, a TM-like 
band of the regular PhC becomes resonant with the cavity mode and 
introduces a new loss channel 38 . We point out that, while Q 

t appears 

to systematically increase with d in the single-mode region, it drops 
rapidly as soon as d increases into the multi-mode region. An ana- 
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Figure 2 | (a)-(b): Electric field (E y ) profiles of the two optimized L3 designs; the holes that were displaced are marked in red. (a): S ix = 0.327a, 
S 2x = 0.257a, S 3st = 0.116a. (b): S lx = 0.337a, S 2x = 0.270a, S 3l = 0.088a, S 4x = 0.323a, S 5x = 0.173a. (c): Dependence of Q f on S 2t and S 3x for S ix = 0.327a; 
the shifts in the design of panel (a) are marked by a white cross, (d): Histograms showing the probability of occurrence of different Q e -values in the 
design of panel (a), for two different disorder magnitudes: a = 0.003a (red) and a = 0.0015a (blue). The black line indicates the ideal Q f . (e): Same as (d), 
for the design of (b). (f): Dependence of Q, on the overall radius and the slab thickness, for the values of S lx — Ss x corresponding to the optimal design in 
(b). On the right of the dashed line, the slab becomes multi-mode at the cavity frequency. 



logous trend of the loss rates as a function of R/a and dla is expected 
for the other cavity designs discussed in this study. In principle, one 
could consider dla and/or R/a as free parameters in the optimization, 
but in that case setting a target wavelength, for a fixed value of d 
that might arise from technological requirements, becomes more 
difficult. 

Often, obtaining the highest possible theoretical Q t is not the main 
goal of optimization. In fact, when Q t gets above a limit set by the 
material and the fabrication process (currently ~ 5 X 10 6 in sil- 
icon 12 ), the experimentally measured Q e is always dominated by 
losses due to disorder and/or absorption, and so weakly affected by 
further increase in Q t . The potential of an automated optimization 
procedure is therefore best exploited when applied to other attractive 
properties. One example consists in maximizing Q t while having the 
smallest possible mode volume, so that Qt/V is as high as possible, 
since the latter is a figure of merit for applications in both cavity 
QED 25,27,28 and non-linear optics 26 . With this in mind, the second 
design we focus on is the HO cavity 6,9 (sometimes also named "zero- 
cell" or "point-shift"), namely the simple defect cavity with the smal- 
lest known mode volume. 

The design of the HO is shown in 3(a); the defining defect is the 
shift of two holes away from each other (S lx ), the thickness of the slab 
taken here is d = 0.5a, while the hole radius is R = 0.25a. For the 
optimization, we also use the consecutive shifts S 2x — S 5x , as well as 
two shifts in the vertical direction, S ly and S 2y . Using Si x , S 2x , S$ x , S ly 
and S 2y , the cavity has already been optimized 6 (following the same 
approach already discussed for the L3 2 ) to a quality factor of Q t = 2.8 
X 10 5 , with a corresponding mode volume V = 0.23(A/«) 3 . Here, we 
improve on this result by on one hand using the genetic optimization, 
and on the other by including S 4x and S 5x . It should be mentioned 
that in the optimization, an allowed range of variation for each para- 
meter is set. For the HO, we find that the maximum allowed S lx is a 



very important parameter, increasing which produces several differ- 
ent optimized designs. All of those are interesting, as increasing S lx 
increases the mode volume but also the Q t of the cavity. More pre- 
cisely, the designs in 3(b)-(d) were obtained by imposing the follow- 
ing restrictions in the genetic algorithm: S lx s 0.25a, S lx s 0.3a, and 
S lx £ OAa, respectively. The ensuing optimal parameters \S\ X , S2 X > 
S 3x , S 4x , S 5x , S ly , S 2y ] are as follows: [0.216a, 0.103a, 0.123a, 0.004a, 
0.194a, -0.017a, 0.067a] (panel (b)); [0.280a, 0.193a, 0.194a, 0.162a, 
0.113a, -0.016a, 0.134a] (panel(c)); and [0.385a, 0.342a, 0.301a, 
0.229a, 0.116a, -0.033a, 0.093a] (panel(d)). The corresponding 
quality factors are Q t = 1.04 X 10 6 (FDTD: 1.04 X 10 6 ), Q t = 1.88 
X 10 6 (FDTD: 1.66 X 10 6 ), and, remarkably, Q t = 8.89 X 10 6 (FDTD: 
8.29 X 10 6 ), while the respective mode volumes are V = 0.25()./n) 3 , V 
= 0.34(a/«) 3 , and V = 0.64(^/«) 3 . The first among these three 
designs (panel (b)) has a mode volume only slightly larger than the 
previous most optimal design 6 , combined to a quality factor almost 
four times larger. The last of the three designs (panel (d)) instead 
shows a more significant increase of the mode volume, but associated 
to an almost 30-fold increase of Q, with respect to the value obtained 
in Ref. 6. The resonance frequencies of the three designs decrease 

caa coa 
with the increase of V and are =0.280, =0.275, and 

2kc 2kc 

=0.269, respectively, while the original cavity with S lx = 0.14 

27TC 

(and no other shifts) of Ref. 9 has =0.292. 

27IC 

Similarly to what we have done above for the L3 cavity, in 
Fig. 3(e)-(g) we present the probability of occurrence of Q e values, 
computed using 1000 random disorder realizations for each design 
and each disorder magnitude. From these histograms it appears 
clearly that, even though design 3 has the highest theoretical Q t /V, 
it might not be the best choice in practice. According to Eq. ( 1 ) in fact, 
depending on the amount of disorder, the maximum value of the 
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Figure 3 | (a): The design of the HO cavity. For quality factor optimization, shifts of the positions of the five neighboring holes in the x-direction and 
the two neighboring holes in the y-direction were introduced, labeled as Si„ S 2x , S 33CI S 4 „ S 5x> S iy and S 2y in the Figure, (b)-(d): Electric field (E y ) 

0.123a, S 4 . 



profiles of three optimized designs with increasing Q and V. (b): S lx = 0.216a, S 2x = 0.103a, S 3 . 
S 2y = 0.067a. (c): S ix = 0.280fl, S 2x = 0.193a, S 3x = 0.194a, S 4x = 0.162a, S 5x = 0.113a, S 



iy ' 



-0.016a, S 



iy 



= 0.004a, S 5x = 
0.134a. (d): Si* 



0.194a, S, 



-0.017a, 



0.385a, S 2x = 0.342a, 



S 3l = 0.301a, S ix = 0.229a, S 5x = 0.116a, Si 7 = —0.033a, S 2y = 0.093a. (e): Histograms showing the probability of occurrence of different Q e -values 
in the design of panel (b), for two different disorder magnitudes: a = 0.003a (red) and a = 0.0015a (blue). The black line indicates the value of Q,. 
(f)-(g): Same as (e), for the designs of (c)-(d). The (7 = 0 line in panel (g) is not visible as Q, occurs beyond the axis boundary. 



actual ratio QJV will in general be achieved for a design having an 
intermediate value of Q t . For example, in the case with a = 0.003a 
(red curves in panels (e)-(g)), the average values of Q e (neglecting 
absorption) computed from the simulations are 3.97 X 10 5 , 5.23 X 
10 5 , and 6.49 X 10 5 , respectively, meaning that the highest QJV 
would in practice be achieved by design 1. On the other hand, for 
the smaller disorder a = 0.0015a (blue curves in panels (e)-(g)), the 
corresponding average values of Q e are 7.22 X 10 5 , 1.12 X 10 6 , and 
2.02 X 10 6 , and the highest average QJ Vis achieved by design 2. For 
both values of a, the expected QJV for the five-shift L3 cavity is lower 
than that for any of the three HO designs, thus the latter should be the 
cavity of choice for applications where QIV is the most important 
figure of merit. We expect that design 3 will dominate for yet smaller 
values of a, which however appear to be currently not achievable in 
practice. Additionally, design 3 could be the preferred choice for 
purely manufacturing reasons if quantum dots are to be inserted in 
the cavity, as the field maxima are further away from the hole edges 
than in the other two cavities. An important final remark is that 
cavities 1 and 2 have already been fabricated 3 *, and experimental 
quality factors consistent with the a = 0.003a disorder histograms 
have been measured (unloaded Q e = 485, 000 for design 2), resulting 
in optical bistability with ultra-low threshold power in silicon. 

Another interesting cavity with potential applications for polar- 
ization-entangled photon generation 43,44 and quantum dot spin read- 
out 45 is the HI 9,13,34 - also known as "single point-defect cavity" - 
formed by one missing hole in the lattice (Fig. 4(a)). The modes of 
this cavity preserve the underlying hexagonal symmetry, and for a 
wide range of parameters, the fundamental (lowest-frequency) res- 
onance is given by two degenerate dipole modes. Based on the electric 
field polarization in the center of the cavity, those are usually referred 
to as the y-polarized (Fig. 4 (b)) and x-polarized (Fig. 4 (c)) modes, 
although in fact the electric field of each of those modes also has a 
non-vanishing component oriented in the orthogonal direction; 
however, since the near-field in the very center of the cavity as well 
as the far-field emission in the z-direction (perpendicular to the slab 



plane) are both truly y (correspondingly x) polarized, this labeling is 
in many cases appropriate, in particular for applications in which a 
quantum dot is placed in the center of the cavity. Here, we optimize 
the quality factor of the dipole mode. We take d = 0.55a and _R = 
0.23a, and the parameters used for design optimization, labeled Si, 
S 2 , and S 3 in Fig. 4(a), are an increase in the side-length of the three 
consecutive hexagonal "rings" around the cavity (which is also equi- 
valent to the increase of the distance from a vertex of a hexagon to the 
cavity center). The previously most-optimal design was achieved 
using only the holes at the vertices of the hexagons (but including 
variations of the hole radii), and has a moderate Q, = 6.2 X 10 4 and a 
mode volume V = 0.47(/ 1 ./n) 38 . Here, using the genetic algorithm 
with the shifts as outlined in Fig. 4(a), we find an optimal design at 
Si = 0.213a, S 2 = 0.070a, and S 3 = -0.009a, with Q t = 1.05 X 10 6 
(FDTD: 0.97 X 10 6 ) and V = 0.62(l/n)\ i.e. we find a 19-fold 
increase in Q t coupled to an increase in V by 32%. For this cavity 
as well, the disorder analysis (panel (d)) suggests that Q e - values close 
to a million can be expected experimentally in state-of-the-art silicon 
systems, i.e. more than an order of magnitude larger than the pre- 
vious theoretical values. The modes lie at a frequency =0.253, 

which is, as expected, slightly lower than that of the unmodified 

a>a 

cavity, =0.270. 

' 2nc 

We note that while the degeneracy of the two dipole modes is an 
attractive feature of the HI, it is lifted by disorder. This is why, in 
panel (e) of Fig. 4, we study the probability of occurrence of the 
splitting between the modes, based on the 1000 disorder realizations 
that were used for the disorder analysis in panel (d). It is important to 
note that there is no absolute way to define an x — y reference frame, 
as three equivalent frames (rotated 60° from one another) exist due to 
the hexagonal symmetry of the cavity. This symmetry is broken if the 
cavity presents preferential orientations of the axes (e.g. introduced 
by lithography). In the case where only random disorder is consid- 
ered, x- and y-polarized modes can turn out to be oriented along 
either of the three x — y reference frames. Thus, what we plot in panel 
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Figure 4 | (a): The design of the HI cavity. The size of the first three hexagonal "rings" of holes is varied for quality factor optimization, with the 
increase of the distance from a hexagon vertex to the center of the cavity given by Si, S 2 and S 3 (marked), (b): Electric field (£„) profile of the y-polarized 
and (c): electric field (E x ) profile of the x-polarized mode, for the optimal design Si = 0.213a, S 2 = 0.070a, S 3 = —0.009a. (d): Histograms showing the 
probability of occurrence of different Q e -values for two different disorder magnitudes: a = 0.003a (red) and a = 0.0015a (blue). The black line 
indicates the ideal Q t . (e): Histograms showing the probability of occurrence of the wavelength splitting AX between the two modes, which are degenerate 
in the disorder-less cavity. 



(e) is only the difference between the resonant wavelengths of the 
higher- and the lower-frequency modes, without any reference to 
their polarization. What is important to note is that the splitting is 
of the order of hundreds of picometers, i.e. two orders of magnitude 
larger than the linewidth corresponding to the typical Q e -s (see top pr- 
axis of panel (d)). This implies that for applications for which overlap 
in frequency between the two modes is needed, some form of post- 
fabrication tuning is required, the possibility for which has already 
been demonstrated 43,46 . This can also be combined with an additional 
modulation of the neighboring holes which increases the emission in 
the vertical direction 47 , which would decrease the Q e and make the 
tuning easier while simultaneously increasing the intensity of the 
collected radiation, which is beneficial for entangled-photon 
generation 44,48 . 

We finally address the hexapole mode of the HI cavity, that typ- 
ically lies at higher frequency than the dipolar modes studied above. 
This mode has previously been optimized to Q t = 1.6 X 10 6 11,13 by 
varying only Si, with S 2 = S 3 = 0 in the sketch of Fig. 4(a). Here, we 
run a global optimization of Q t by varying the three shifts, with d = 
0.55a and _R = 0.22a. We could improve the previous value to 
Q t = 3.2 X 10 6 (FDTD: 3.1 X 10 6 ), obtained for the optimal values 

Si = 0.271a, S 2 = 0.039a, and S 3 = 0.018a. The resonance frequency 
coa 

of this mode is =0.261, while in the unmodified cavity of the 

2kc ' 
same d and R, this mode is not present. 

The figures of merit of the seven designs that were optimized 
above are summarized in Table I. Notice that for the higher disorder 
a i, the values of (Q e ) do not differ dramatically among the different 
designs. This reflects the fact that, as has been shown in a previous 
analysis 37 , in the disorder- dominated regime (when Qd^>Qt), the 
losses are nearly design-independent. 

Discussion 

The designs obtained here for the three most widespread PhC defect- 
cavities consistently show that the quality factor of these cavities can 
be systematically optimized to well above 10 s by adjusting only a few 
structural parameters (only shifts of hole positions were used here), 



with small increases in the mode volumes (within 50% with the 
exception of the third HO design) as compared to those of the cor- 
responding non-optimized designs. In carrying out our analysis, we 
have tried to include the radii of holes next to the cavity as additional 
free parameters, but this brought no significant improvement. We 
therefore restricted to shifts of hole positions only, as these are easily 
controlled in the fabrication process. Our scheme leaves the possibil- 
ity open to use the hole radii as free parameters for independent 
optimization of some additional figure of merit. 

A very important conclusion of our analysis is that the parameter 
space of such structural variations has to be explored globally, using 
an automated optimization tool. We find that the genetic algorithm 
is an excellent tool to handle this task, even when seven parameters 
are included in the computation, as was the case for the HO cavities. 
Employing this algorithm was possible only due to the computa- 
tional advantage of the GME - to compute the number of cavity 
configurations that was needed for the optimization using a first- 
principle tool like FDTD or a Finite-Element Method (FEM) would 
require either an enormous computational power, or time of the 
order of years. This computational advantage made it also possible, 
for each cavity type, to vary more structural parameters than in 
previous optimization works, thus bringing to an even larger increase 
in the quality factors. Finally, we note that the GME returns not only 
the quality factor but also the full mode profile of the cavity modes, 
thus the same procedure which was presented here can be used for 
optimization of different quantities depending on the practical 
requirements, like Q/V for applications in non-linear optics, or 
optimal electric field profiles for e.g. sensing technologies 23 , which 
we expect to explore in the next future. 

The statistical analysis of Q-factors including structural disorder 
shows clearly that the designs obtained here are as robust to disorder 
as other ultrahigh- Q designs 37 . In particular, for the L3 cavity with a 
theoretical Q-factor Q, = 5.1 X 10 6 , state-of-the-art fabrication qual- 
ity in silicon should easily result in experimental Q-factors around 2 
X 10 6 , in the same range as current ultrahigh- Q designs based on a 
PhC waveguide 1012 . On the other hand, our experience shows that 
systematic structural variations can rapidly suppress the Q-factor of 
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Table I | Summary of the simulated resonant frequency, quality factor, mode volume, and Q/Vratio for the seven optimal designs studied 
in this work. The last four columns list average values computed over 1 000 disorder realizations, in which both hole radii and position 
are subject to Gaussian fluctuations with standard deviation a\ = 0.003a and 02 = 0.001 5a 
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an optimal structure. Hence, an important conclusion of our analysis 
is that, whenever a design needs to be significantly varied (e.g. when 
strongly modifying the ratios R/a and/ or d/a, or the refractive index n 
in order to, e.g., operate at a different wavelength), a new optimiza- 
tion must be carried out in order to obtain the best structural design 
adapted to the new requirements. 

The range of applicability of the present scheme is very broad. We 
are currently applying it to several systems, including the optimiza- 
tion of the Q-factor of PhC cavities with low index contrast 49 , of 
cavities based on a thicker slab (d/a > 1, typically needed to host 
some quantum nanostructures 50 ), but also to the maximization of the 
trapping capabilities of slot cavities designed for biological sensing 23 , 
and to the simultaneous optimization of first and second order 
modes of a cavity built in a PhC structure with doubly resonant 
bandgap for efficient second-harmonic generation 51 . In all these 
cases, the global optimization brought to very promising results. A 
further application might include the optimization of PhC structures 
based on a silicon-on-insulator design 52,53 . Perhaps even more 
importantly, our work has general implications about the future 
research in photonic crystals, as it shows that the domain of possible 
structural designs is still largely unexplored, while simultaneously 
demonstrating how its exploration can be done efficiently. 

Methods 

The guided-mode expansion. The GME method used in this work 38 is based on 
expanding the mode of a cavity on the basis of the guided modes of an effective 
dielectric slab. In the z-direction (orthogonal to the slab plane), the whole space from 
minus to plus infinity is included analytically; in the slab plane, a finite super-cell is 
taken and periodic boundary conditions are assumed. The super-cell size in our 
simulations ranges from L x = 16a to L x = 20a in the ^-direction and from 

L y = 12(v / 3a/2) to Ly — 20 f \/3a 1 2^ in the y-direction, depending on the cavity. 

The supercell naturally defines a grid of in-plane wave-vectors G = 2n(n x /L x , n y IL y ) 
(with integer n x and n y ), each of them labeling a wave guided inside the slab and 
evanescent along the z-direction. The expansion is carried out onto a truncated set of 
N G such waves (centered at G = 0) and the number N G is increased until convergence 
is reached (typically JV G > 5000). For all the optimizations, only the fundamental (a = 
1 in Ref. 38) guided mode is used, since for the thickness assumed here, the admixture 
of higher guided modes is negligible at the cavity resonant frequencies. In the multi- 
mode region of Fig. 2(f), the second (y. — 2) guided mode is also included. To correct 
for a small systematic offset in the quality factor Q t due to the periodic boundary 
conditions, an averaging of the losses of the Bloch modes of the cavity over a grid in k- 
space containing 10 to 20 points is performed. With these numerical parameters, 
obtaining the Q for a given cavity configuration takes of the order of five minutes 
using one processor and 2 GB of RAM. When the slab becomes multi-mode, the 
distribution of losses in fc-space presents a sharp, narrow peak around the k values for 
which the cavity mode is exactly in resonance with the TM-like band of the regular 
PhC 38 . In that case a finer grid {or numerical integration) in /c-space is needed to 
correctly compute the quality factor, and simulating one cavity configuration takes of 
the order of a few hours with the same computational resources. 

FDTD simulations. For the FDTD simulations, we use the freely available MIT 
Electromagnetic Equation Propagation (MEEP) package 40 . Typical simulations 

parameters are supercell size 40a in the ;c- direction and 30^^^ / 2j in the y- 

direction, and spatial resolution 20/(3. Taking advantage of the parallelized (MPI) 
version of the software, a converged computation takes of the order of ten hours on a 
16-core processor with 32 GB of RAM. A more detailed comparison of 



computational time and accuracy between GME, FDTD and the Finite-Element 
Method (FEM) can be found in Ref. 37. 

Genetic optimization. The implementation of the genetic algorithm included in the 
Global Optimization Toolbox (MATLAB and Global Optimization Toolbox Release 
2012b, The Math Works, Inc., Natick, Massachusetts, United States), which starts 
from a random initial population (i.e. a set of points in parameter space) and goes on 
to create a sequence of generations (new sets of such points) where the 'fittest 5 
individuals are kept. The algorithm incorporates an array of evolutionary inspired 
techniques, including cross-over, random mutations, and natural selection of 
individuals. An 'individual' in our case is simply one Q-computation for a particular 
set of cavity parameters, and the higher the quality factor, the higher the assigned 
'fitness' of the individual (i.e. its probability to survive to the next generation and/or to 
be mixed with another individual to produce an 'offspring' lying in parameter space 
somewhere in between the two). With the increase of the number of free parameters, 
the number of generations needed for convergence increases. In the presented work, 
the maximum is reached for the HO cavity with 7 parameters, for which —300 
generations each consisting of 120 individuals are needed for convergence. This can 
however be greatly improved if a rough optimization is first carried out (with a large 
allowed range for the free parameters), followed by a finer optimization centered 
around the rough maximum. The longest optimization we ran thus took about a week 
on twelve CPUs with 32 GB of RAM. This would have taken more than ten years to 
finish using the same machine but employing an FDTD solver instead of the GME. 

Mode volume. Here we adopt the definition of the mode volume that is most 
commonly used to quantify the performance of a nanocavity: 



This definition is the one needed in cavity quantum electrodynamics as a measure of 
the radiation -matter coupling with a pointlike two-level system placed at the position 
where the field intensity has a maximum. Other definitions are more suited to 
different purposes and might in some cases differ dramatically from the one above 11 . 
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